Project2

Gaussian random field with application of INLA

Part I Multivariate normal distribution

Let \(\boldsymbol{x} = (x_1, \dots, x_n), n = 100\) be multivariate normal distributed with \(E(x_i) = 0, Var(x_i) = 1\), and \(Corr(x_i, x_j) = e^{-0.1|i - j|}\)

1. Compute and image the covariance matrix \(\boldsymbol{\Sigma}\) of \(\boldsymbol{x}\)

2. Find the lower Cholesky factor \(\boldsymbol{L}\), such that \(\boldsymbol{L}\boldsymbol{L}^T = \boldsymbol{\Sigma}\), of this covariance matrix, and image.

3. Sample \(\boldsymbol{x} = \boldsymbol{Lz}\), where \(\boldsymbol{z}\) is a length n random vector of independent standard normal variables. Plot the sample.

4. Find the precision matrix \(\boldsymbol{Q}\) of the covariance matrix, and compute the lower Cholesky factor \(\boldsymbol{L}_Q\), such that \(\boldsymbol{L}_Q\boldsymbol{L}_Q^T = \boldsymbol{Q}\), of this matrix. Image these matrices and compare them to the images obtained in a) and b)

5. Sample \(\boldsymbol{x}\) by solving \(\boldsymbol{L}_Q^T\boldsymbol{x} = \boldsymbol{z}\), where \(\boldsymbol{z}\) is a length n random vector of independen standard normal variables. Plot the sample.

6. Permute the ordering of variables in \(\boldsymbol{x}\), and redo the exercises.

---

Solution to Part I

a)

Given that \(\boldsymbol{\Sigma} = e^{-0.1|i - j|}\). The covariance matrix can be expressed as follows: \[\begin{equation*} \Sigma = \begin{pmatrix} 1 & e^{- 0.1 h_{12}} &\dots & e^{- 0.1 h_{1n}} \\ e^{- 0.1 h_{21}} & 1 & \dots & e^{- 0.1 h_{2n}} \\ \vdots & \vdots & \ddots & \vdots \\ e^{- 0.1 h_{n1}} & e^{- 0.1 h_{n2}} & \dots & 1 \end{pmatrix} \end{equation*}\]

### b) According to the cholesky decomposition rule, \(\boldsymbol{L}\) is the lower triangular matrix for \(\boldsymbol{\Sigma}\), it can be easily computed from R using L = chol(Sigma). It is then plotted as below.

c)

Sample using \(\boldsymbol{x} = \boldsymbol{L}\boldsymbol{z}\) transforms the zero-mean, standard normal random variales to the random variables with the desired covariance matrix.

d)

The precision matrix \(\boldsymbol{Q}\) is the inverse of the covariance matrix \(\boldsymbol{\Sigma}\), it is computed using Q = solve(Sigma) in R. The three matrices are thereby depicted as follows. Since the covariance matrix is not singular, given that the it belongs to the Matern family, thus it is analytically guaranteed to have positive definite property. Therefore, both precision matrix and the lower triangular precision matrix exist.

e)

Similarly, the expected random samples can be generated using the inversion of the above formula, thus \(\boldsymbol{L}_Q^T\boldsymbol{x} = \boldsymbol{z}\)

f)

Permute \(\boldsymbol{x}\) to make randomise the ordering of the grid, the asscociated covariance matrix can be thereby modified in a sparse way.

Part II Gaussian random fields and Kriging

The purpose of this computer exercise is to give an introduction to parameter estima- tion and kriging for Gaussian random field models for spatial data.

We assume the following observation model on the unit square:

\[\begin{equation*} y(\boldsymbol{s}_j) = x(\boldsymbol{s}_j) + \epsilon_j, \ \ \ j = 1, \dots, N, \end{equation*}\]

where \(\epsilon_j ~ \mathop{N}(0, \tau^2)\) are independent measurement noise terms. Further, consider a Matérn covariance function for the Gaussian random field \(\boldsymbol{x(s)}\):

\[\begin{equation*} Cov(x(\boldsymbol{s}_i), x(\boldsymbol{s}_j)) = \Sigma_{i, j} = \sigma^2(1 + \phi h)\exp(-\phi h), \end{equation*}\]

where \(h\) denotes the Euclidean distance between the two sites \(\boldsymbol{s}_i\) and \(\boldsymbol{s}_j\).

We assume the mean increases with east and north coordinates as follows: \(\mu_j = \alpha((s_{j1} - 0.5) + (s_{j2} - 0.5))\), for site \(\boldsymbol{s}_j = (s_{j1}, s_{j2})\) on the unit square.

---

2.1 Simulation

Simulate \(N = 200\) random sites in the unit square and plot them. Form the covariance matrix using \(\sigma = 1, \phi = 10, \tau = 0.05\). Take its Cholesky decomposition and simulate dependent zero-mean Gaussian data variables, then add the mean using \(\alpha = 1\). Plot your observations.

The true mean of the field is expressed as \[\begin{equation*} \mu_i = \alpha((s_{i1} - 0.5) + (s_{i2} - 0.5)) = -\alpha + \alpha s_{i1} + \alpha s_{i2} = \underbrace{ \begin{bmatrix} 1 & s_{i1} & s_{i2} \\ \end{bmatrix}}_{\boldsymbol{h}^T(\boldsymbol{s}_i)} \underbrace{ \begin{bmatrix} -\alpha \\ \alpha \\ \alpha \end{bmatrix}}_{\boldsymbol{\beta}} \end{equation*}\] where \(s_{i1}, s_{i2}\) are the location from east and north direction in the grid.

---

2.2 Paramter estimation

We will now use the simulated data to estimate the model parameters \(\alpha, \sigma^2, \tau^2, \phi\) using maximum likelihood estimation. Iterate between the update for the mean parameter, and updating the covariance parameters. Monitor the likelihood function at each step of the algorithm to check convergence.

The mean of the field is modelled by \(p(\boldsymbol{x})\) and the imperfect information \(\boldsymbol{y} = (y_1, \dots, y_m)\) conditional on \(\boldsymbol{x}\) can be modelled by \(p(\boldsymbol{y}|\boldsymbol{x})\), which can be expressed as follows: \[\begin{equation*} p(\boldsymbol{x}) = \mathop{N}(\boldsymbol{H}\beta, \boldsymbol{\Sigma}), \quad p(\boldsymbol{y}|\boldsymbol{x}) = \mathop{N}(\boldsymbol{Fx}, \boldsymbol{T}) \end{equation*}\] Therefore, the marginal likelihood of the data is \[\begin{equation*} p(\boldsymbol{y}) = \mathop{N}(\boldsymbol{G\beta}, \boldsymbol{C}), \quad \boldsymbol{G} = \boldsymbol{FH}, \quad \boldsymbol{C} = \boldsymbol{F\Sigma F^T} + \boldsymbol{T} \end{equation*}\] The log-likelihood as a function of \(\beta\) and unknown fixed nuisance parameters \(\boldsymbol{\theta}\) in the prior covariance matrix \(\boldsymbol{\Sigma} = \boldsymbol{\Sigma(\boldsymbol{\theta})}\), and/or the likelihood noise matrix \(\boldsymbol{T} = \boldsymbol{T}(\boldsymbol{\theta})\) becomes \[\begin{equation*} l(\boldsymbol{\theta}, \boldsymbol{\beta}) = -\frac{m}{2} \log(2\pi) - \frac{1}{2}\log|\boldsymbol{C}| - \frac{1}{2}(\boldsymbol{y} - \boldsymbol{G}\boldsymbol{\beta})^T\boldsymbol{C}^{-1}(\boldsymbol{y} - \boldsymbol{G\beta}) \end{equation*}\]

The MLEs of \(\boldsymbol{\beta}\) and \(\boldsymbol{\theta}\) are obtained by \[\begin{equation*} (\hat{\boldsymbol{\beta}}, \hat{\boldsymbol{\theta}}) = \underset{\boldsymbol{\beta}, \boldsymbol{\theta}}{\arg\max}\{l(\boldsymbol{\beta, \theta})\} \end{equation*}\] For fixed \(\boldsymbol{\theta}\), the MLEs of \(\boldsymbol{\beta}\) can be determined analytically. \[\begin{equation*} \frac{dl}{d\boldsymbol{\beta}} = \boldsymbol{G}^T\boldsymbol{C}^{-1}\boldsymbol{y} - \boldsymbol{G}^T\boldsymbol{C}^{-1}\boldsymbol{G}\boldsymbol{\beta} = \boldsymbol{0}, \quad \hat{\boldsymbol{\beta}} = (\boldsymbol{G}^T\boldsymbol{C}^{-1}\boldsymbol{G})^{-1}\boldsymbol{G}^T\boldsymbol{C}^{-1}\boldsymbol{y} \end{equation*}\] Whereas for fixed \(\boldsymbol{\beta}\), the MLE of nuisance parameters \(\boldsymbol{\theta}\) can be obtained by numerical maximization. Let \(\boldsymbol{z} = \boldsymbol{y} - \boldsymbol{G\beta}\), and \(\boldsymbol{Q} = \boldsymbol{C}^{-1}\). For each component of \(\boldsymbol{\theta}_r\), \(r = 1, \dots, d\), in this case, \(\theta\) has 3 components (\(\sigma^2, \eta, \tau^2\)). The score of the log-likelihood becomes \[\begin{equation*} \frac{dl}{d\theta_r} = -\frac{1}{2}\text{trace}(\boldsymbol{Q}\frac{d\boldsymbol{C}}{d\theta_r}) + \frac{1}{2}\boldsymbol{z}^T\boldsymbol{Q}\frac{d\boldsymbol{C}}{d\theta_r}\boldsymbol{Qz} \end{equation*}\]

The above mentioned score can be solved iteratively using Fisher scoring algorithm. To achieve the numerical stability, the expected Hessian is applied, which is \[\begin{equation*} \mathop{E}(\frac{d^2l}{d\theta_rd\theta_{\bar{r}}}) = -\frac{1}{2}\text{trace}(\boldsymbol{Q}\frac{d\boldsymbol{C}}{d\theta_{\bar{r}}}\boldsymbol{Q}\frac{d\boldsymbol{C}}{d\theta_r}) \end{equation*}\]

## [1] "1.16100274094671  , iter no is  0"
## [1] "0.0983478023616246  , iter no is  1"
## [1] "0.0637476756371519  , iter no is  2"
## [1] "0.0330323876752982  , iter no is  3"
## [1] "0.0216765191888463  , iter no is  4"

## 
## Estimated sigma is  0.785 ; True sigma is  1 
## Estimated phi is  12.72 ; True phi is  10 
## Estimated tau is  -0.06224 ; True tau is  0.05 
## Estimated beta1 is  -1.47 ; True beta1 is  -1 
## Estimated beta2 is  1.63 ; True beta2 is  1 
## Estimated beta3 is  1.88 ; True beta3 is  1NULL

---

2.3 Kriging

We will now use the estimated model parameters to perform kriging prediction. Predict variables \(x(s)\), where predictions sites lie on a regular grid of size 25x25 for the unit square. Visualize the Kriging surface and the prediction standard error. Compare with the true field.

## [1] "The prediction error is  0.303326130373502"

---

Part III Integrated nested Laplace Approximations (INLA)

3.1 Simple Linear Regression

In the last part of this exercise we explore the R-INLA package.

The ski jumping data set contains 26 observations of measured lengths in ski jumping competitions (in meters) between 1961 and 2011.

Visualisation of the ski jumping data set

In Figure we observe a clear (almost linear) trend in the measured jumping lengths to increase with the years.

Therefore, we use linear approach for the statistical modelling of this data, where the years \(x_i\) are the covariates and the lengths \(y_i\) are the responses for \(i=1,\dots,26\): \[\begin{equation*} \mathbb{E}[y_i] = \mu + \beta x_i,\quad \mathbb{Var}[y_i]=\tau^{-1}. \end{equation*}\]

This can be posed as a latent Gaussian model suited for the INLA framework.

1. The response depends on the linear predictors \(\eta\) as \(y|x,\theta = \Pi \pi(y_i|\eta_i,\tau)\) where \(\pi(y_i|\eta_i,\tau)\sim\mathcal{N}(\eta_i,\tau^{-2})\)
2. The parameter of the linear predictor \(eta_i=\mu+x_i\beta\) are independent Gaussian with some variance and mean zero.
3. The hyperparameter \(\tau\) is only one-dimensional.

## 
## Call:
##    "inla(formula = Length ~ Year, data = skiData)" 
## Time used:
##     Pre = 4.73, Running = 0.494, Post = 0.259, Total = 5.48 
## Fixed effects:
##                  mean      sd 0.025quant  0.5quant 0.975quant      mode kld
## (Intercept) -4029.646 106.615  -4240.649 -4029.650  -3818.934 -4029.646   0
## Year            2.126   0.054      2.019     2.126      2.232     2.126   0
## 
## Model hyperparameters:
##                                          mean   sd 0.025quant 0.5quant
## Precision for the Gaussian observations 0.072 0.02      0.038     0.07
##                                         0.975quant  mode
## Precision for the Gaussian observations      0.116 0.066
## 
## Expected number of effective parameters(stdev): 2.00(0.00)
## Number of equivalent replicates : 13.00 
## 
## Marginal log-Likelihood:  -89.58

Posterior marginals for selected effects

## Mean            3.83993 
## Stdev           0.552976 
## Quantile  0.025 2.93482 
## Quantile  0.25  3.44617 
## Quantile  0.5   3.77761 
## Quantile  0.75  4.16357 
## Quantile  0.975 5.09958

FB: What does all this mean?…..

3.2 GLMM with random effects

Last, we use INLA to analyse the “Seeds” data set.

## 
## Call:
##    c("inla(formula = formula, family = \"binomial\", data = data, Ntrials 
##    = n, ", " control.predictor = list(compute = TRUE), control.family = 
##    list(link = \"logit\"))" ) 
## Time used:
##     Pre = 5.06, Running = 0.456, Post = 0.251, Total = 5.77 
## Fixed effects:
##               mean    sd 0.025quant 0.5quant 0.975quant   mode kld
## (Intercept) -0.429 0.115     -0.656   -0.429     -0.204 -0.428   0
## x1          -0.272 0.156     -0.580   -0.272      0.033 -0.271   0
## x2           1.066 0.146      0.782    1.066      1.353  1.065   0
## 
## Random effects:
##   Name     Model
##     plate IID model
## 
## Model hyperparameters:
##                         mean       sd 0.025quant 0.5quant 0.975quant mode
## Precision for plate 18930.69 20071.40      30.63 12701.76   73765.97 7.83
## 
## Expected number of effective parameters(stdev): 3.25(1.16)
## Number of equivalent replicates : 6.46 
## 
## Marginal log-Likelihood:  -73.75 
## Posterior marginals for the linear predictor and
##  the fitted values are computed